20/04/21
For manuscript notes & figures.
Based on http://127.0.0.1:8888/lab/tree/Scratch/XeF2-notebooks/xe-xef2_comp_120221-dist.ipynb and SO stuff from http://127.0.0.1:8888/lab/tree/Scratch/XeF2-notebooks/XeF2_ePS-expt_comp_121120_4d-dev.ipynb
CURRENT VERSION (21/04/21): copied in some Lyx LaTex, set cell metadata and tested. Working... not very slick. Need to rerun for PDF figures.
12/02/21
Comparison of Xe & XeF2 data + experiments. No SO treatment herein.
Basic conclusions: atomic Xe parameters (cross-sections and $\beta$) look broadly similar to XeF2 results, except:
For XeF2 only plus SO treatment, see previous notebook XeF2_ePS-expt_comp_271020_4d_v111120-dist.html.
!hostname
!conda env list
import sys
import os
from pathlib import Path
import numpy as np
# import epsproc as ep
import xarray as xr
import matplotlib.pyplot as plt
from datetime import datetime as dt
timeString = dt.now()
# For module testing, include path to module here, otherwise use global installation
if sys.platform == "win32":
modPath = r'D:\code\github\ePSproc' # Win test machine
winFlag = True
else:
modPath = r'/home/femtolab/github/ePSproc/' # Linux test machine
winFlag = False
sys.path.append(modPath)
import epsproc as ep
# Plotters
from epsproc.plot import hvPlotters
# Multijob class dev code
from epsproc.classes.multiJob import ePSmultiJob
hvPlotters.setPlotters(width = 1000, snsStyle="whitegrid")
# import bokeh
# import holoviews as hv
# hv.extension('bokeh')
# For class, above settings don't take, not sure why, something to do with namespaces/calling sequence?
# Overriding snsStyle does work however... although NOT CONSISTENTLY????
# AH, looks like ordering matters - set_style LAST (.set seems to override)
import seaborn as sns
sns.set(rc={'figure.figsize':(10,6)}) # Set figure size in inches
sns.set_context("paper")
sns.set_style("whitegrid") # Set plot style
sns.set_palette("Paired") # Set colour mapping
# sns.set_theme(context='paper', style='whitegrid') # Not valid for v0.9.0
# Try direct fig type setting for PDF output figs
from IPython.display import set_matplotlib_formats
# set_matplotlib_formats('png', 'pdf')
set_matplotlib_formats('svg', 'pdf')
# # Scan for subdirs, based on existing routine in getFiles()
# fileBase = r'D:\projects\ePolyScat\xef2\XeF2_highRes_wf' # Data dir on Stimpy
# fileBase = r'D:\VMs\Share\ePSshare\xe\Xe_wf' # Data dir on Stimpy
# jobDirs = [r'D:\VMs\Share\ePSshare\xe\Xe_wf\orb19_4d'] # Check OK for single dir specified method
# Set explicit list of dirs
jobDirs = [r'D:\VMs\Share\ePSshare\xe\Xe_wf\orb19_4d', r'D:\projects\ePolyScat\xef2\XeF2_highRes_wf\orb21_A1G',
r'D:\projects\ePolyScat\xef2\XeF2_highRes_wf\orb22_E1G', r'D:\projects\ePolyScat\xef2\XeF2_highRes_wf\orb24_E2G']
data = ePSmultiJob(jobDirs = jobDirs, verbose = 0)
# keys = [4,5,6] # Set for 4d data only
# data.scanFiles(keys = keys)
data.scanFiles() # All files
data.jobsSummary()
(Note: for the atomic case, there is an HU continuum symmetry allowed, but the matrix elements are, effectively, zero - this triggers the warning above.)
# Change labels for plotting
# data.data['orb19_4d']['jobNotes']['orbLabel'] = 'Xe atomic 4d, S/Ag'
for key in data.data.keys():
data.data[key]['jobNotes']['orbLabel'] = data.data[key]['jobNotes']['orbLabel'].split(',')[-1]
# print(data.data[key]['jobNotes']['orbLabel'].split(',')[-1])
# Set branching ratios
# Stack XS data to new data structure
# NOTE: this is not quite correct, since it forces all data to same Ehv axis, but OK for a quick hack, and easy to pull out branching ratios
dsXS = xr.Dataset() # Set blank dataset, this is easier for stacking, probably
lText = []
for key in data.data.keys():
if key is not 'orb19_4d': # Skip Xe atomic case
subset = data.data[key]['XS'].sel({'Sym':'All', 'XC':'SIGMA'}) # Set XS data, all syms only
# NOTE currently missing full dataset resolution for orb24, so try interp. (2eV not 1eV step size)
# Note dropna to ensure no NaNs, see http://xarray.pydata.org/en/latest/interpolation.html#interpolating-arrays-with-nan
if key == 'orb24_E2G':
subset = subset.dropna('Eke').interp(Eke = dsXS.Eke, method = 'cubic')
# dsXS[key] = subset.copy() # Set .copy() for safety here!
dsXS[data.data[key]['jobNotes']['orbLabel']] = subset.copy() # Set .copy() for safety here!
# Convert to Xarray & normalise
dsXS = dsXS.to_array().rename({'variable':'channel'}).squeeze()
dsXS['total'] = dsXS.sum('channel') # Sum over channels
# Normalise...
dsXS = dsXS/dsXS['total']
# Load experimental data
dataPathExpt = Path(r'D:\projects\XeF2_Soleil_2019\RF_data_analysis_021020')
# Set data attribs in dict similar to ePS results structure
dataFiles = {}
dataFiles['SIGMA'] = {'Data':'XS', 'File':r'XeF2_Xe4d_4d32_4d52_cross_sec_all_photon_energies_02102020.dat',
'States':['$4d_{3/2}$', '$4d_{5/2}$']}
dataFiles['BETA'] = {'Data':'Beta', 'File':r'XeF2_Xe4d_beta_all_photon_energies_02102020.dat',
'States':['$\Pi_{1/2} (4d_{3/2})$', '$\Delta_{3/2} (4d_{3/2})$', '$\Sigma_{1/2} (4d_{5/2})$',
'$\Pi_{3/2} (4d_{5/2})$', '$\Delta_{5/2} (4d_{5/2})$']}
dataFiles['BR-All'] = {'Data':'Branching ratios', 'File':r'XeF2_Xe4d_branching_all_photon_energies_all_states_01112020.dat',
'States':['$\Pi_{1/2} (4d_{3/2})$', '$\Delta_{3/2} (4d_{3/2})$', '$\Sigma_{1/2} (4d_{5/2})$',
'$\Pi_{3/2} (4d_{5/2})$', '$\Delta_{5/2} (4d_{5/2})$']}
dataFiles['BR-SO-summed'] = {'Data':'Branching ratios SO summed', 'File':r'XeF2_Xe4d_branching_all_photon_energies_SO_av_01112020.dat',
'States':['$\Pi$', '$\Delta$', '$\Sigma$']}
# Update with ion data
dataFiles2 = {}
dataFiles2['ION-low'] = {'Data':'XS', 'File':r'XeF2_ion_yield_low_energy_cal.txt'}
dataFiles2['ION-high'] = {'Data':'XS', 'File':r'XeF2_ion_yield_high_energy_cal.txt'}
# Update with branching ratios
# dataFilesBR = {}
# dataFilesBR['All'] = {'Data':'Branching ratios', 'File':r'XeF2_Xe4d_branching_all_photon_energies_all_states_01112020.dat',
# 'States':['\pi1/2 (4d3/2)', '\delta3/2 (4d3/2)', '\sigma1/2 (4d5/2)', '\pi3/2 (4d5/2)', '\delta5/2 (4d5/2)']}
# dataFiles['SO-summed'] = {'Data':'Branching ratios SO summed', 'File':r'XeF2_Xe4d_branching_all_photon_energies_SO_av_01112020.dat',
# 'States':['\pi', '\delta', '\sigma']}
# Read data files and convert to Xarray
# 27/10/20 added quick hack to set 2nd array for ion yield data
dataList = []
dataList2 = []
for key in dataFiles:
dataIn = np.loadtxt(dataPathExpt/dataFiles[key]['File'])
# Convert to Xarray
dataXr = xr.DataArray(dataIn[:,1:], dims=('Ehv','State'), coords={'Ehv':dataIn[:,0], 'State':dataFiles[key]['States'][0:dataIn.shape[1]-1]})
dataXr.attrs = dataFiles[key]
dataXr.attrs['dataPath'] = dataPathExpt
dataXr.name = key
dataList.append(dataXr)
# Stack to Xarray
dataExpt = xr.concat(dataList, "XC")
dataExpt['XC'] = list(dataFiles.keys())
for key in dataFiles2:
dataIn = np.loadtxt(dataPathExpt/dataFiles2[key]['File'])
# Convert to Xarray
dataXr = xr.DataArray(dataIn[:,1].squeeze(), dims=('Ehv'), coords={'Ehv':dataIn[:,0]}) # Only 1D datasets in this case
dataXr.attrs = dataFiles2[key]
dataXr.attrs['dataPath'] = dataPathExpt
dataXr.name = key
dataList2.append(dataXr)
dataExpt2 = xr.concat(dataList2, "XC")
dataExpt2['XC'] = list(dataFiles2.keys())
Here set Seaborn pallete for colour-mapping, and include labels etc. - aim to be publication ready.
# sns.color_palette("cubehelix")
# sns.set_palette("Reds")
# sns.set_palette("Paired")
# Compare with computational results
marker = None #'x'
Etype = 'Ehv'
Erange = [70, 255]
pType = 'SIGMA'
pltObj, lTextComp = data.plotGetCroComp(pType=pType, Etype = Etype, Erange = Erange, returnHandles = True)
lText = lTextComp.copy()
# Add expt data - cross-secions
scale = dataExpt.sel({'XC':pType}).max() / data.dataSets['orb19_4d']['XS'].sel({'XC':pType, 'Type':'L'}).max() # Set scaling to match computational data
(dataExpt.sel({'XC':pType})/scale).dropna('State').plot.line(x='Ehv', marker=marker, ls='--');
# Add expt data - ion yields
ionData = 'ION-high'
scale = dataExpt2.sel({'XC':ionData}).max() / data.dataSets['orb21_A1G']['XS'].sel({'XC':pType, 'Type':'L'}).max() # Set scaling to match computational data
(dataExpt2.sel({'XC':ionData})/scale).plot.line(x='Ehv', marker=marker, ls='--');
# Manual legend fix
lText.extend(dataExpt.sel({'XC':pType}).dropna('State').coords['State'].data)
lText.extend([ionData])
plt.legend(lText)
# Fix axis labels
if pType == 'SIGMA':
plt.ylabel('Cross-section/Mb')
# plt.title('Cross-sections (normalised to computational results)')
plt.title('')
else:
plt.ylabel(r"$\beta_{LM}$")
# plt.title(r"$\beta_{LM}$")
# Fix plot x-axis
plt.xlim(Erange);
# plt.yscale("log") # Use log scale y-axis?
Comparison of experimental (dashed lines) and computational (solid lines) cross-sections. Computational cross-sections are in Mb, and experimental results are scaled to match the peak value.
# Beta comparison plot over orbs
pType = 'BETA'
pltObj, lText = data.plotGetCroComp(pType=pType, Etype = Etype, Erange = Erange, returnHandles = True)
# Add expt data
dataExpt.sel({'XC':pType}).dropna('State').plot.line(x='Ehv', marker=marker, ls='--');
# Manual legend fix
lText.extend(dataExpt.sel({'XC':pType}).dropna('State').coords['State'].data)
plt.legend(lText, loc='upper right')
# Fix axis labels
if pType == 'SIGMA':
plt.ylabel('XS/Mb')
plt.title('XS (normalised to computational results)')
else:
plt.ylabel(r"$\beta_{LM}$")
# plt.title(r"$\beta_{LM}$")
plt.title('')
# Fix plot x-axis
plt.xlim(Erange);
Comparison of experimental (dashed lines) and computational (solid lines) $\beta$ parameters.
# Beta comparison plot over orbs
Erange = [120, 200]
pType = 'BETA'
pltObj, lText = data.plotGetCroComp(pType=pType, Etype = Etype, Erange = Erange, returnHandles = True)
# Add expt data
dataExpt.sel({'XC':pType}).dropna('State').plot.line(x='Ehv', marker=marker, ls='--');
# Manual legend fix
lText.extend(dataExpt.sel({'XC':pType}).dropna('State').coords['State'].data)
plt.legend(lText, loc='upper right')
# Fix axis labels
if pType == 'SIGMA':
plt.ylabel('XS/Mb')
plt.title('XS (normalised to computational results)')
else:
plt.ylabel(r"$\beta_{LM}$")
plt.title(r"$\beta_{LM}$")
# Fix plot x-axis
plt.xlim(Erange);
# Compare with computational results
Erange = [75, 250]
pType = 'BR-SO-summed'
marker = 'x'
# pltObj, lText = data.plotGetCroComp(pType=pType, Etype = Etype, Erange = Erange, returnHandles = True)
dsXS.swap_dims({'Eke':'Ehv'}).sel(**{Etype:slice(Erange[0], Erange[1])}, Type='L').plot.line(x='Ehv');
# Add expt data - cross-secions
# scale = dataExpt.sel({'XC':pType}).max() / data.dataSets['orb21']['XS'].sel({'XC':pType, 'Type':'L'}).max() # Set scaling to match computational data
scale = 1
(dataExpt.sel({'XC':pType})/scale).dropna('State').plot.line(x='Ehv', marker=marker, ls=':');
# Manual legend fix
lText = lTextComp.copy()[1:]
lText.extend(dataExpt.sel({'XC':pType}).dropna('State').coords['State'].data);
plt.legend(lText, loc='upper right');
plt.ylabel("Branching ratio")
plt.title("");
Comparison of experimental (dashed lines) and computational (solid lines) branching ratios. Experimental results are summed over spin-orbit bands to allow direct comparison with the computational results.
From http://127.0.0.1:8888/lab/tree/Scratch/XeF2-notebooks/XeF2_ePS-expt_comp_121120_4d-dev.ipynb
Just final bits here, hopefully.
from epsproc.geomFunc.geomUtils import genllLList
from epsproc.geomFunc.geomCalc import w3jTable
# Gen QNs for specific (L,S) case
L = 2
S = 0.5
Llist = np.array([[L,L+S,S], [L, L, S], [L,L-S,S]]) # Note this needs to be 2D array in current form of function
QNs = genllLList(Llist, uniqueFlag = False)
# Then calc 3js....
backend = 'sympy'
form = 'xdaLM' # '2d' # 'xdaLM' # xds
nonzeroFlag = True
dlist = ['L', 'J', 'S', 'Lambda', 'Omega', 'Sigma'] # Set dims for reference
thrj = w3jTable(QNs = QNs, nonzeroFlag = nonzeroFlag, form = form, dlist = dlist, backend = backend)
# And primed terms (will be identical at this point, but set dims for multiplication later)
thrjP = w3jTable(QNs = QNs, nonzeroFlag = nonzeroFlag, form = form, dlist = [item + 'p' for item in dlist], backend = backend)
# Reformate by comsolidating +/- terms as modulation for XS and square
thrjUS = thrj.unstack('mSet').fillna(0)
# thrjXSmod = 2*(thrjUS.sum('Sigma').where((thrjUS['Lambda']>-0.5)&(thrjUS['Omega']<0.5), drop=True))**2 # Note Lambda & Omega anti-phase!
thrjXSmod = 2*(thrjUS.sum('Sigma').where((thrjUS['Lambda']>-0.5)&(thrjUS['Omega']<0.5), drop=True))**2 # Note Lambda & Omega anti-phase!
thrjXSmod['Omega'] = np.abs(thrjXSmod['Omega'])
pdTable, _ = ep.multiDimXrToPD(thrjXSmod, colDims = 'Lambda')
pdTable
# SO branching ratios
# 11/11/20 - converted to function
def simulateSOBR(dsXS, thrjSOTerm, BRtype = 'All', coupling = 'Lambda', stateLabels = ['$\Sigma$', '$\Pi$', '$\Delta$']):
# Set Lambda terms
# dsXS['Lambda'] = dsXS['channel']
# dsXS['Lambda'] = [0.0, 1.0, 2.0]
dsXSO = dsXS.assign_coords({coupling:('channel',thrjSOTerm[coupling])}).swap_dims({'channel':coupling}) # Assign Lambda for multiplication
dsXSO = dsXSO.assign_coords({f'{coupling}-Labels':(coupling, np.asarray(stateLabels)[dsXSO[coupling].astype(int).data])}) # Assign labels, based on int values as indexers (CRUDE!)
dsXSO = dsXSO * thrjSOTerm
# dsXS
# Convert to branching ratios (total)
if BRtype == 'All':
# dsXSO['total'] = dsXSO.sum(['Lambda', 'Omega', 'lSet']) # Sum over channels
dsXSO['total'] = dsXSO.sum(thrjSOTerm.dims)
elif BRtype == 'J':
dsXSO['total'] = dsXSO.sum(['Lambda', 'Omega']) # Sum over channels, except J - GIVES same result?
elif BRtype == 'None':
dsXSO['total'] = 1.0
else:
print(f'BRtype={BRtype} not recognised.')
return
# # Set branching ratios
dsXSO = dsXSO/dsXSO['total']
return dsXSO
Label levels by ($\Omega$, $\Lambda$), with no additional summation over terms. Q: is this justified/realistic?
pType = 'BR-All'
# Selected term(s) only - LOOKS GOOD FOR coupling='Omega' case, but opposite phases for coupling='Lambda'
dsXSO = simulateSOBR(dsXS, thrjXSmod, BRtype = 'None', coupling='Omega') # This looks good
# dsXSO = simulateSOBR(dsXS, np.abs(thrjXSmodCoherent), BRtype = 'None', coupling='Lambda') # Testing coherent case
# Unnorm
# dsXSO *= dsXSO['total']
# Selected term(s) only
J=1.5
OLset = [[2,1.5], [1,0.5]] # [Omega, Lambda] pairs to match expt.
# OLset = [[1,0.5],[2,1.5]] # [Omega, Lambda] pairs to match expt.
# OLset = [[1.5,2], [0.5,1]] # Check ordering OK!
# J=2.5
# OLset = [[2,2.5], [1,1.5], [0,0.5]] # [Omega, Lambda] pairs to match expt.
# subset = xr.DataArray()
subset = []
for OL in OLset:
subset.append((dsXSO).swap_dims({'Eke':'Ehv'}).sel(**{Etype:slice(Erange[0], Erange[1])}, Type='L', J=J, Lambda=OL[0], Omega=OL[1]))
subsetXS = xr.concat(subset, dim='Omega').squeeze()
# Renorm
subsetXS = subsetXS/subsetXS.sum('Omega')
subsetXS.plot.line(x=Etype)
# Selected states only
if J == 1.5:
statesPlot = ['$\\Delta_{3/2} (4d_{3/2})$', '$\\Pi_{1/2} (4d_{3/2})$']
else:
statesPlot = ['$\\Delta_{5/2} (4d_{5/2})$', '$\\Pi_{3/2} (4d_{5/2})$', '$\\Sigma_{1/2} (4d_{5/2})$']
(dataExpt.sel({'XC':pType, 'State':statesPlot})/scale).plot.line(x='Ehv', marker=marker, ls=':');
# Manual legend fix
# lText = list(subsetXS['Omega'].data)
# lText = list(OLset)
lText = [f"{item.split(' ')[0]}$(sim)" for item in statesPlot]
# lText = list(subsetXS['Lambda'].data)
lText.extend(statesPlot)
# lText.extend(dataExpt.sel({'XC':pType}).dropna('State').coords['State'].data)
plt.legend(lText, loc='lower right')
plt.ylabel("Branching ratio")
# plt.title(f'Simulated branching ratios, J={J}, Hund\'s case (c), with subselection (no summation)')
plt.title("");
plt.ylim([0.3, 0.68]); # Set ylims to avoid large peaks
Comparison of experimental (dashed lines) and computational (solid lines) branching ratios for $J=3/2$.
# Selected term(s) only
# Selected term(s) only - LOOKS GOOD FOR coupling='Omega' case, but opposite phases for coupling='Lambda'
# For J = 5/2 'Lambda' case is possibly better, at least for state ordering? DO NOT UNDERSTAND THIS YET!
# scale = 1
dsXSO = simulateSOBR(dsXS, thrjXSmod, BRtype = 'None', coupling='Omega')
# Unnorm
# dsXSO *= dsXSO['total']
# Selected term(s) only
# J=1.5
# OLset = [[2,1.5], [1,0.5]] # [Omega, Lambda] pairs to match expt.
J=2.5
# OLset = [[2,2.5], [1,1.5], [0,0.5]] # [Omega, Lambda] pairs to match expt.
OLset = [[0,0.5], [1,1.5], [2,2.5]] # [Omega, Lambda] pairs to match expt - REVERSED STATE ORDERING, justified by ion state energy ordering?
# subset = xr.DataArray()
subset = []
for OL in OLset:
subset.append((dsXSO).swap_dims({'Eke':'Ehv'}).sel(**{Etype:slice(Erange[0], Erange[1])}, Type='L', J=J, Lambda=OL[0], Omega=OL[1]))
subsetXS = xr.concat(subset, dim='Omega').squeeze()
# Renorm
subsetXS = subsetXS/subsetXS.sum('Omega')
subsetXS.plot.line(x=Etype)
# Selected states only
if J == 1.5:
statesPlot = ['$\\Delta_{3/2} (4d_{3/2})$', '$\\Pi_{1/2} (4d_{3/2})$']
else:
statesPlot = ['$\\Delta_{5/2} (4d_{5/2})$', '$\\Pi_{3/2} (4d_{5/2})$', '$\\Sigma_{1/2} (4d_{5/2})$']
(dataExpt.sel({'XC':pType, 'State':statesPlot})/scale).plot.line(x='Ehv', marker=marker, ls=':');
# Manual legend fix
# lText = list(subsetXS['Omega'].data)
# lText = list(OLset)
lText = [f"{item.split(' ')[0]}$(sim)" for item in statesPlot]
# lText = list(subsetXS['Lambda'].data)
lText.extend(statesPlot)
# lText.extend(dataExpt.sel({'XC':pType}).dropna('State').coords['State'].data)
# plt.legend(lText, loc='upper right')
plt.legend(lText, loc='lower right')
plt.ylabel("Branching ratio")
# plt.title(f'Simulated branching ratios, J={J}, Hund\'s case (c), with subselection (no summation)')
plt.title('');
plt.ylim([0.1, 0.5]); # Set ylims to avoid large peaks
Comparison of experimental (dashed lines) and computational (solid lines) branching ratios for $J=5/2$.
Interactive^ plots per orbital & symmetry. In these plots the solid lines show the 'mixed' gauge results, and dashed lines show length and velocity gauges as bounds on the values.
^ Use controls in top right of plot. Also data sets can be turned on/off using the legend entries.
Note that Xe atomic case is the first plot below, ignore the out-of-bounds $\beta$s for (null) HU symmetry case, these are artefacts. Also of note in this case is that the individual symmetry channels show very little variation with energy, hence most of the overall (symmetry-summed) structure is form intra-channel effects (interferences).
# XC & Betas
data.plotGetCro(pType='SIGMA', Etype = Etype, Erange = Erange, backend = 'hv')
As previously for XeF2 case; for atomic Xe there are 5 degenerate channels (4d components) labelled by 'it'.
data.lmPlot(dataType = 'matE', Erange = Erange, Etype = Etype, thres = 0.1, logFlag = True, selDims = {'Type':'L', 'it':1}, fillna = True, cmap = "vlag", figsize = (15,8))
# cmap = "cubehelix" ) # Quite good for showing resonance features.
data.lmPlot(dataType = 'matE', Erange = Erange, Etype = Etype, thres = 0.1, logFlag = False, selDims = {'Type':'L', 'it':1}, fillna = True, cmap = "vlag", pType = "phase", figsize = (15,8))
# cmap = "cubehelix" ) # Quite good for showing resonance features.
import scooby
scooby.Report(additional=['epsproc', 'xarray', 'jupyter'])